#include "refraction.h"

#include "vector3.h"
#include <cmath>
#include <cstdlib>

namespace {
  const uint32_t kRussianRouletteDepth = 3;
}

namespace raytrace {

Refraction::Refraction(float refract_index,
                       Vector3Arg reflectance_color,
                       Vector3Arg emittance_color)
  : Material(reflectance_color, emittance_color),
    refract_index_(refract_index) {
}
  
void Refraction::GetOutDirection(Vector3Arg in, Vector3Arg normal,
                                 uint32_t depth,
                                 std::vector<float>& ipdfs,
                                 std::vector<Vector3>& outs) const {
  float ndi = Vector3::dot(normal, in);
	Vector3 reflect_dir = Vector3::sub(in, Vector3::scale(normal, 2 * ndi));
	reflect_dir.normalize();
	
	const Vector3 nl = Vector3::scale(normal, (ndi < 0 ? 1.f : -1.f));
	
	const bool into = (Vector3::dot(normal, nl) > 0);
	const float nc = 1.0f; // refractive index of air
	const float nt = refract_index_;
	const float nnt= into ? nc/nt : nt/nc;
	const float ddn = Vector3::dot(in, nl);
	const float cos2t = 1 - (nnt * nnt * (1 - ddn * ddn)); 
	if (cos2t < 0) {
		// Total internal reflection
    outs.push_back(reflect_dir);
    ipdfs.push_back(1.0f);
    return;
	}

  Vector3 transmit_dir = Vector3::scale(in, nnt);
  transmit_dir = Vector3::sub(transmit_dir, 
                              Vector3::scale(nl, (ddn * nnt + sqrtf(cos2t))));
  transmit_dir.normalize(); 
  
  const float a = nt - nc;
  const float b = nt + nc;
  const float R0 = (a * a) / (b * b);
  const float c = (into ? -ddn : Vector3::dot(transmit_dir, normal));
  const float c1 = 1 - c; 
  
  const float Re = R0 + ((1 - R0) * c1 * c1 * c1 * c1 * c1);
  const float Tr = 1 - Re;
  const float P = 0.25f + 0.5f * Re;
			
  if (depth > kRussianRouletteDepth) {
    // Russian roulette
    if(drand48() < P) {
        outs.push_back(reflect_dir);
        ipdfs.push_back(Re / P);
    } else {
        outs.push_back(transmit_dir);
        ipdfs.push_back(Tr / (1.0f - P));
    }

    return;
  } 
//		float thetaIncident = acosf(ddn);
//		float thetaTransmission = acosf(c);
//		const float fresnelTransmit =
//          Fresnel::CalculateDielectric(nnt, thetaIncident, thetaTransmission);
  const float fresnelTransmit = 1.0;
  
  outs.push_back(reflect_dir);
  ipdfs.push_back(2 * M_PI * Re);
  
  outs.push_back(transmit_dir);
  ipdfs.push_back(2 * M_PI * Tr * fresnelTransmit);
}                                      

}
